Recovering boundary conditions in inverse Sturm-Liouville problems 

Norbert Rohrl 

Abstract. We introduce a variational algorithm, which solves the classical inverse Sturm-Liouville 
problem when two spectra are given. In contrast to other approaches, it recovers the potential as well as 
the boundary conditions without a priori knowledge of the mean of the potential. Numerical examples 
show that the algorithm works quite reliable, even in the presence of noise. A proof of the absence 
of strict local minimizers of the functional supports the observation, that a good initial guess is not 
essential. 



1. Introduction 

The inverse Sturm-Liouville problem was first systematically studied by Borg in 1946 l t . He already 
proved that all information needed to reconstruct the potential is in two sequences of eigenvalues, and 
applied this to the question if one could hear the mass density of a guitar string. 

With modern computers, the question for efficient algorithms to actually compute the potential 
gained importance [7j- This work was inspired by two different approaches to this problem. 

One was developed by Rundell and Sacks |10j and is based on the Gelfand-Levitan-Marchenko kernels, 
ft is elegant and efficient, but also invariably needs the mean J Q Qdx of the potential and the boundary 
conditions as additional inputs besides the two spectra. The second method is variational and was created 
by Brown, Samko, Knowles, and Marietta |2j. It does not need the mean as a separate input, but it is 
unknown if it also can be used to recover the boundary conditions. 

We also want to mention a related recovery method by Lowe, Pilant, and Rundell which uses a 
finite basis ansatz, and solves the inverse problem by Newton's method without requiring the mean as 
input. 

In this paper we extend the variational method we introduced in [2] to recover potential and boundary 
conditions in the case when only two finite sequences of eigenvalues are given. Having less reliable 
information it is not as robust under noisy input, but we still get reasonable results. 

In the following section we will define the functional and exhibit some essential properties. The 
numerical examples are discussed in the third section and section four finally contains the proof of the 
absence of strict local minimizers. 

2. Definition and Properties of the Functional 

We consider the Sturm-Liouville equation 

(SL) — u" + q(x)u = \u 

on [0, 1] with q(x) G L 2 ([0, 1],K) real, and separated boundary conditions 

{h hi) h u(0) + u'(0) = 0, hiu(l) + u'(l) = 0. 
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The corresponding eigenvalues satisfy the asymptotic formula [2] 



(2.3) 



A, 



7r 2 n 2 + 2(hi — h ) + / q(s)ds + a n , 



where (o„) € I 2 . It is a classical result Q15|, that the potential q is uniquely determined by two sequences 
of eigenvalues corresponding to boundary conditions (hohi) and (/10/12) with h\ ^ h 2 . Moreover it can 
be shown, that those sequences also uniquely determine the boundary conditions py. 

Therefore, two sequences of eigenvalues contain all information necessary to recover the potential as 
well as the corresponding boundary conditions. For notational convenience we write the parameters of 
the two Sturm-Liouville problems as vectors 



(ho,hi,q), q 2 ■= {h ,h 2 ,q), 



Aq .n for the n-th eigenvalue of problem q i: and 

q := (h 0) hi,h 2) q) 

for the full problem. 

Now we define a least squares functional on the eigenvalues, which has the solution of the inverse 
problem as zero. 

Definition 2.1. Suppose we are given (partial) spectral data Xq.^ with (i,n) in 7 C {1, 2} x N of 
an unknown Sturm-Liouville problem Q = (77q, 77i, 772, Q)- For a trial problem q and positive weights 



U); 



we define the functional 



(2.4) 



(i,n)el 



We note that G(q) is positive, and zero if and only if both given sequences of eigenvalues match those 
of q. If we have full knowledge of the two sequences Xq^u, (h n) € {1, 2} x N, this determines q uniquely, 
and hence q = Q. 

To find such a q, we minimize the functional with a conjugate gradient descent algorithm. First, for 
numerical stability it is good to know, that for each pair of interlacing sequences, 



Al,„ < \ 2 ,n < Ai jW +i 



or 



Aa,n < Al, n < A2, n +1 i 



which satisfy the asymptotics (|2.3|l . there is a q, with G(q) = [3], 
The gradient of X qi , n wrt. q { is |3] 



/ 9\ qi ,„ \ 

1 dh > 



VA n 



dh 2 
\ dq J 



( 



9q L ,n 



(1)^,1 



(.2 



<n(l)4 
V 9q i<n ( x ) J 



= 1. 



where g qi ,n denotes the eigenfunction corresponding to \ qi , n with \\g q 
It follows that the gradient of the functional is given by 

VG(q) = 2 w i,n(Ag 4 ,„ - A Q . in )VAq., n , 

(i,n)e/ 

if (u>i t n) is summable. 

Theorem 14. 21 below shows that the gradients VAq.,„ are linearly independent in R 3 x 7 2 (0, 1). This 
immediately implies the essential convexity of the functional: 
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Theorem 2.2. Let q l = (ho, hi, q) and q 2 = (ho, hi, q) be two Sturm- Liouville problems with h\ 7^/12- 
/// is finite or (tOi^ n ) is summable, the functional G has no local minima at q with G(q) > 0, i.e. 

VG(g) = G(q) = . 

Thus a conjugate gradient algorithm will not get trapped in local minima, as we will also observe in 
the examples. 

3. Numerical Examples 

We use the standard Polak-Ribiere conjugate gradient descent algorithm |5] to approximate the 
gradient flow and thus minimize the functional. To give the basic idea, we explain the simpler steepest 
descent: 

(i) choose initial potential and boundary conditions q^ = (h^\ h^\ h^\ q^) 

(ii) while G(qV>) too big do 

(a) compute the gradient VG(q^) 

(b) minimize the one dimensional function G(q^) — aVG(q^) wrt. a 

(c) set equal to the minimizing potential 

This straight forward minimization scheme has a major disadvantage: consecutive gradients are always 
orthogonal. To avoid this, conjugate gradient descent computes the direction for the one dimensional 
minimization using the current and previous gradients. 

Note that the boundary points <?q.,n(0) and g q . in (l), needed in the computation of the gradient 
VAq..„, can be computed in a numerically well behaved way. Given the eigenvalue A q . in , we can compute 
a multiple of the eigenfunction by solving an initial value problem. The value of the eigenfunction at 
the boundary then is just the (exactly known) initial value divided by the L 2 norm of the initial value 
solution. 

We first apply the algorithm to a popular non-continuous potential |2j E] 

Q(x) = (7x - 0.7)x(o.i,o.3](z) + (3-5 - 7x)x (0 .3,o.5] (x) + 4x(o.7,o.9]( a + ^X(x.o.9,i](x) , 
and choose 

(3.1) Q=(3,3,0,Q) and g (0) = (2, 4, -1, 0) 

with / = {1, 2} x {0, . . . , 29} and weights uii >n = 1. This will be our default setting, unless noted otherwise. 
In the figures we write A2 = ||<7 — QH2 for the L 2 error of the reconstruction and be = (ho; hi; hi) for the 
boundary conditions of the current approximation. 

The results (figure ^) are comparable to the alternative boundary condition example with given 
boundary conditions in [2]. We get quite good results at 150 iterations, which keep getting better as we 
minimize the functional. 

From around 150 iterations the boundary conditions stay almost constant. This suggests to reset 
the trial potential q^ to zero at some iteration number j, while keeping the boundary conditions. In 
other words, we reset the forth component of = (h^ , q^) , and keep the others fixed. In 

figure we indeed attain a significantly better approximation by setting the potential to zero a couple 
of times. In all our examples this was a very useful strategy to get faster convergence. But it is just 
heuristics - we do not really know how to choose the optimal number of iterations j. In practice, we wait 
until the boundary conditions stabilize and then set the potential to zero. This can be repeated until the 
convergence speed of the functional G does not improve any more. 

The graph in figure demonstrates, that the reconstruction of a smooth potential, using the same 
boundary conditions and number of eigenvalues, yields more accurate results. 
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Figure 1. Graph of q^ 150 \ q( 840 > (light) and boundary conditions, G(q), and L 2 error 
versus the number of iterations. 
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Figure 2. Above example with setting g( 30 ) = 
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,(60) 



q 



(no) = o. 



This overall behavior is also true for worse guesses of the initial boundary conditions. If we take for 
example 

= (0,0,0,0), 

the boundary conditions converge slowly, but steadily (figure0J. Again, setting qV> = a couple of times 
increases the speed of convergence dramatically. 
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Figure 3. Approximation of a smooth potential. 




Figure 4. The convergence of the boundary conditions and of the functional for initial 
problem = (0,0,0,0). 



In the case of noisy data, the algorithm is of course much more unstable than the version with fixed 
boundary conditions [S]. For the following examples, we add random noise |AQ. jra — Xq.,u\ < r to the 
eigenvalues. To see the limitations of this approach, we first use r = 0.1 and Hq = ho = Hi = hi = 3, 
H 2 = hi = (figure [5j, i.e. we already start with the correct boundary conditions. 

The algorithm passes through a potential, which is reasonably close to the original potential Q. But 
from there, the steepest descent leads to a potential, which is not in any way similar to the one we want to 
recover. In the graph of the L 2 error, we see that there are roughly 20 iterations of good approximations. 
Afterwards the approximations quickly get worse than our initial guess. 

Yet, for smaller errors in the eigenvalues, this effect is less dramatic. Setting r = 0.01 and using the 
problem i|3.1|) . we get good approximations for around 100 iterations (figure EJ. The L 2 error also rises 
more slowly than for the case with larger noise level. 

Setting «?( 89 ) = again improves the performance significantly (figure 01. We get reasonable approx- 
imations for all 400 iterations but the first 20 after each setting the potential to zero. 
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Figure 5. Example with errors of magnitude r = 0.1. The best L 2 approximation at 
23 iterations and q( 55 \ 



4. Linear Independence of the Gradients 

First, we borrow the central lemma of the independence proof of the original functional [5]. To that 
end, we define the Wronskian [/, g] = fg' — f'g and the bilinear form 

r : i^QO, 1],M) 2 — > E 

(f,9) ^ f [f,9]dx 

which is bounded by 

ir^mi/MMU, i.e. ||r(/, .)|| = 11/11^. 



(We use the definition ||/||ffi = y ||/||f,2 + 11/11^2 with distributional derivatives.) 

Let Si,n,(j and Ci^ n ^ q be the solutions of the differential equation (|SL|) for the eigenvalue parameter 
Aq., n and initial values 

= -hi, <„, g (l) = -hi- 

Lemma 4.1 Given two Sturm- Liouville problems q 1 = (ho,hi,q) and q 2 — (ho,Ii2,q) with 

hi ^ hi, we have 

T{ci,n,qSi, n ,q,g 2 q . tm ) = {-l) l (h 2 - hi)5 n , m 8i^ for all i,j G {1,2} and m,n e N 
for the normalized eigen] 'unctions g q ., n and Si >n ^ q , Ci^ n ^ q as defined above. 
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Using the alternative bilinear form 

f: i? 1 ([0, 1], K) x (M 3 x £ 2 ([0, 1],K)) — > M 

(f,(a,b,c,g)) ~ -2^f'gdx + f(l)b + f(l)c + f(0)a ' 

and integration by parts 

T(f,g) = -2 / 1 /'.9dx + /.g(l)-/ 5 (0), 
Jo 

we get the corresponding statement for the gradients X q ^ m 

^{Ci,n,qSi,n,q:^ Xq^ra) — ^{^i,n,q^i,n,qj 9q^ ,m) ( ^) (^2 ^l)^n,m^i,j • 

Finally, since f is bounded by 

|f (/, (a,b,c,g))\ < 2\\f\\ Hl \\g\\ L 2 + V^II/IMM + l&l + |c|) < 2\\f\\ Hl (\\g\\ L3 + \a\ + \b\ + |c|) 
and, in particular, continuous in the 2nd component, we immediately get the linear independence theorem. 
Theorem 4.2. With the notations of the above lemma, the set of gradients of the eigenvalues 

{V\, n |(i 1 n)6{l,2}xN} 

is linearly independent in M 3 x L 2 . 
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Figure 7. Example with errors of magnitude r — 0.01 and qr' 89 ) set to zero. The best 
L 2 approximation at 237 iterations and g( 389 ) . 
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Proof. Suppose for some fixed (i, n) we have 



VAq.^„ — ^ QfcVAfc 



in M 3 x L 2 , where e R and VAfc = VA„ m with {jk,mk) ^ {i,n). But this would imply 

(-l) l (/l 2 - hi) = f (Cj^^Sj^,,, VA,.. n ) = f Ci^ q Si^ q , ^ a k^ J = ^ f (c i;ni9 S i;ni9 , ttfeVAfc) = 

V feGN / feGN 
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